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We present a multiscale simulation algorithm for amorphous materials, which we illustrate and 
validate in a canonical case of dense granular flow. Our algorithm is based on the recently proposed 
Spot Model, where particles in a dense random packing undergo chain-like collective displacements 
in response to diffusing "spots" of influence, carrying a slight excess of interstitial free volume. We 
reconstruct the microscopic dynamics of particles from the "coarse grained" dynamics of spots by 
introducing a localized particle relaxation step after each spot-induced block displacement, simply 
to enforce packing constraints with a (fairly arbitrary) soft-core repulsion. To test the model, 
we study to what extent it can describe the dynamics of up to 135,000 frictional, viscoelastic 
spheres in granular drainage simulated by the discrete-element method (DEM). With only five 
fitting parameters (the radius, volume, diffusivity, drift velocity, and injection rate of spots), we find 
that the spot simulations are able to largely reproduce not only the mean flow and diffusion, but also 
some subtle statistics of the flowing packings, such as spatial velocity correlations and many-body 
structural correlations. The spot simulations run over 100 times faster than DEM and demonstrate 
the possibility of multiscale modeling for amorphous materials, whenever a suitable model can be 
devised for the coarse-grained spot dynamics. 



I. INTRODUCTION 

The geometry of static sphere packings is an age-old 
problem [l| with current work focusing on jammed ran- 
dom packings 0, , but how do random packings flow? 
Here, we consider the case of granular drainage |^ , which 
is of practical importance (e.g. in pebble-bed nuclear re- 
actors |3, and also raises fundamental questions in 
non-equilibrium statistical mechanics Q. In fast, dilute 
flows, Boltzmann's kinetic theory of gases can be modi- 
fied to account for inelastic collisions p, but slow, dense 
flows (as in Fig. require a different description due 
to long-lasting, many-body contacts |9|- Although bal- 
listic motion may occur at the nano-scale 10] (< 0.01% 
of a grain diameter), collisions do not result in random 
recoils, as in a gas. 

In crystals, diffusion and flow are mediated by de- 
fects, such as vacancies and dislocations, but in disor- 
dered phases it is not clear what, if any, "defects" might 
facilitate structural rearrangements. Perhaps the only 
candidate in the literature is an empty "void" in the 
random packing into which a single particle may hop, 
thereby displacing the void. The void mechanism was 
proposed by Eyring for viscous flow Ull and has re- 
appeared in theories of the glass transition [l^ , shear flow 
in metallic glasses [l3|, compaction in vibrated granular 
materials Q, and granular drainage from a silo [TEj . 
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(b) 

FIG. 1: (Color online) A simulation of the experiment in 
Ref. p by discrete element simulations, (a) First, 55,000 
glass beads are poured into a quasi-two-dimensional silo (8 
beads deep) and let come to rest, (b) Slow drainage occurs 
after a slit orifice is opened. (The grains are identical, but 
colored by their initial height.) 



but it is now seen as unrealistic. In glasses, coopera- 
tive relaxation (involving; many particles at once) has 
been observed |T?ill7|. presumably facilitated by free vol- 
ume [H| l2Cj|. In granular drainage, the Void Model 
gives a reasonable fit to the mean flow 0, H3] , and yet 
it grossly over-predicts diffusion |£j . 

A collective mechanism for random-packing dynamics 
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(a) (b) (c) 

FIG. 2: (Color online) The mechanism for structural rearrangement in the Spot Model. The random displacement r s of a 
diffusing spot of free volume (dashed circle) causes affected particles to move as a block by an amount r p (a), followed by an 
internal relaxation with soft-core repulsion (b), which yields the net cooperative motion (c). (The displacements, typically 100 
times smaller than the grain diameter, are exaggerated for clarity.) 



has recently been proposed to resolve this paradox and 
applied to granular drainage [23|. The basic hypothesis, 
shown in Fig. Efa), is that a block of neighboring grains 
makes a small, correlated downward displacement, 

Ar p = -wAr s , (1) 

in response to the random upward displacement, Ar s , of 
a diffusing "spot" of free volume. The coefficient w (more 
generally, a smooth function of the particle-spot separa- 
tion) is set by local volume conservation. In the simplest 
approximation, a spot carries a slight excess of interstitial 
volume, V s , spread uniformly across a sphere of radius 
R s . When the spot engulfs N particles, each of volume 
V p , the model predicts w ~ V 3 /NV P ~ Acfi/cb 2 , where Ac/) 
is the local change in volume fraction cf>. Allowing for 
some spot overlaps yields the estimate w « 10 -2 — 10 -3 
from the observation that Acp/cj) w 1% in dense flows, 
which is consistent with diffusion measurements in exper- 
iments H,!^ and our simulations below. Unlike the Void 
Model (which requires w = 1), each grain's "cage" of 
nearest neighbors also persists over long distances |9| ; the 
Spot Model is able to capture such features of drainage 
experiments, while remaining simple enough for mathe- 
matical analysis, because it does not explicitly enforce 
packing constraints, only the tendency of nearby parti- 
cles to diffuse together. 

In order to preserve valid packings, a multiscale spot 
algorithm has also been suggested 23], which we imple- 
ment here for the first time. As shown in Fig. each 
spot-induced block displacement (a) is followed by a re- 
laxation step (b), in which the affected particles and their 
nearest neighbors experience a soft-core repulsion (with 
all other particles held fixed). The net displacement in 
(c) involves a cooperative local deformation, whose mean 
is roughly the block motion in (a). It is not clear a priori 
that this procedure can produce realistic flowing pack- 
ings, and, if so, whether the relaxation step dominates 
the simple dynamics from the original model. 

To answer these questions, we calibrate and test the 
Spot Model against large-scale computer simulations of 



granular drainage, shown in Fig. [I] Simulations are ad- 
vantageous in this case since three-dimensional packing 
dynamics cannot easily be observed experimentally. We 
begin by running discrete-element method (DEM) simu- 
lations, described in section^ We then calibrate the free 
parameters in the Spot Model by measuring various sta- 
tistical quantities from the DEM simulation, as described 
in IIIII In section HVl we describe the computational im- 
plementation of the Spot Model, before carrying out a 
detailed comparison to DEM in section IVl 



II. DEM SIMULATION METHOD 

We employ a DEM plEsj to simulate N frictional, 
visco-elastic, spherical glass beads of diameter, d = 3mm, 
mass m under the influence of gravity g = 9.81ms -1 . 
Similar to the experiments of Refs. [y, [23 the silo has 
width bOd and thickness 8d with side walls at x = ±25d 
and front and back walls at y = ±4<i, all with friction 
coefficient \i = 0.5. The initial packing is generated by 
pouring N = 55,000 particles in from a fixed height of 
z = 170d and allowing them to come to rest under grav- 
ity, filling the silo up to H Q w 110<i. We also studied a 
taller system with N = 135,000 generated by pouring 
particles in from a height of z = 495<i, which fills the silo 
to H Q w 230<i. We refer to these systems by their initial 
height H Q . Drainage is initiated by opening a circular 
orifice of width Sd centered at x = y = in the base of 
the silo (z = 0). A snapshot of all particle positions is 
recorded every 2 x 10 4 time steps (St = 1.75 x 10 _6 s). 
Once particles drop below z = — lOd, they are removed 
from the simulation. 

The particles interact according to Hertzian, history 
dependent contact forces. If a particle and its neighbor 
are separated by a distance r, and they in compression, 
so that S = d — |r| > 0, then they experience a force F = 
F n + Ft, where the normal and tangential components 
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are given by 



F„ = 



IJfd (-kt As f 
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Here, n = r/|r|. v n and are the normal and tan- 
gential components of the relative surface velocity, and 
k n j and 7 n? £ are the elastic and viscoelastic constants, 
respectively. As*; is the elastic tangential displacement 
between spheres, obtained by integrating tangential rel- 
ative velocities during elastic deformation for the life- 
time of the contact, and is truncated as necessary to 
satisfy a local Coulomb yield criterion |F^| < /z|F n |. 
Particle- wall interactions are treated identically, but the 
particle- wall friction coefficient is set independently. For 
the current simulations we set k t = |fc n , and choose 
k n = 2 x 10 5 mg/d. While this is significantly less than 
would be realistic for glass spheres, where we expect 
k n ~ 10 10 mg/d, such a spring constant would be pro- 
hibitively computationally expensive, as the time step 

— 1/2 

must have the form St oc k n for collisions to be mod- 
eled effectively. Previous simulations have shown that 
increasing k n does not significantly alter physical results 
|2q . We make use of a time step of St = 1.75 x 10 _6 s, 
and damping coefficients 7 n = 7t = 



III. CALIBRATION OF THE MODEL 

We first look for evidence of spots in the DEM simu- 
lation and then proceed to calibrate the model. All the 
calibrations are carried out for the small H Q = llOd sys- 
tem, after which the same parameters are used for the 
larger H Q = 230d system. 

The theory predicts large numbers of spots (since many 
are released as each particle exits the silo), so we seek a 
statistical signature of the passage of many spots. We 
therefore consider the spatial correlation for velocities in 
the x direction, defined by 



C(r) 



(u x (0)u x (r)) 
V(MO) 2 ) (u x (r) 2 ) 



where the expectations are taken over all pairs of veloc- 
ities (u x (0),u x (r)) of particles separated by a distance r 
in a given test region. For a uniform spot influence out to 
a cutoff radius, R s > d, as shown in Fig. Efa), two ran- 
dom particle displacements are either identical, if they 
are caused by the same spot, or independent. In that 
case, the spatial velocity correlation function is given by 



C(r) 



3 r 
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r < 2R S 







h (jk) ' - (4) 
r>2R s 



which is the intersection volume of spheres of radius R s 
separated by r (scaled to 1 at r = 0). The shape of C(r) 
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FIG. 3: (Color online) Comparison of velocity correlations 
calculated over the time period 0.52s < t < 1.57s. Calcu- 
lations are based on particle velocity fluctuations about the 
mean flow in a 16d x 16d region high in the center of the 
container. For H Q — IIOgL 



is affected by the relaxation step in Fig. EJb), but the 
decay length is set by the spot size. 

As shown in Fig. 03 we see spatial velocity correlations 
in the DEM simulations at the scale of several particle 
diameters, consistent with the spot hypothesis. Similar 
correlations have also been seen in experiments 26] using 
the methods of Choi et al. 9], which attests to the gen- 
erality of the phenomenon, as well as the realism of the 
simulations. Since the shape of C(r) is not precisely that 
of Eq. (@J, due to relaxation effects, we fit the simulation 
data to a simple decay, C(r) = ae~ r ^ with (3 = 1.87d. 
We also fit a simple decay of the same form to Eq. Q. 
finding f3 = 0.72i? s , so we infer R s — 2.60<i as the spot ra- 
dius. Thus a grain has significant dynamical correlations 
with neighbors up to three diameters away. 

Next, we infer the dynamics of spots, postulating inde- 
pendent random walks as a first approximation. We as- 
sume that spots drift upward at a constant mean speed, 
v s = Az s /At : (determined below), opposite to gravity, 
while undergoing random horizontal displacements of size 
Ax s in each time step At. The spot diffusion length, 
b s = Var(Ax s )/2Az s , is obtained from the spreading of 
the mean flow away from the orifice. In DEM simula- 
tions, the horizontal profile of the vertical velocity com- 
ponent is well described by a Gaussian, whose variance 
grows linearly with height, as shown in Fig.0] Applying 
linear regression gives Var(i^) = 2.28zd+ 1.60<i 2 , which 
implies b s = 2.28d/2 = 1.14d. To reproduce the spot 
diffusion length, we chose Az s = O.ld and Ax s = 0.68<i. 

The typical excess volume carried by a spot can now 
be obtained from a single bulk diffusion measurement. 
From Eq. JIJ, the particle diffusion length, b p: is given 

by 
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FIG. 4: (Color online) Comparison of the mean velocity pro- 
file, for three different heights calculated over the time period 
4.37s <t< 5.25s once steady flow has been established. The 
Spot Model successfully predicts a Gaussian velocity profile 
near the orifice and the initial spreading of the flow region 
with increasing height, although the DEM flow becomes more 
plug- like higher in the silo. 
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We measure b p in the DEM simulation by tracking the 
variance of the x displacements of particles that start 
high in the silo as a function of their distance dropped. 
We find b p = 2.86 x 10~ 3 d and thus w = 2.50 x 10~ 3 . Dur- 
ing steady flow in the DEM simulation, a typical pack- 
ing fraction of particles is 57.9%, so a spot with radius 
R s = 2.60d influences on average 81.7 other particles. 
Thus we find that a spot carries roughly 20% of a parti- 
cle volume: V s = 81.7V p /w = 0.205V P . 

The three spot parameters so far (radius, R s , diffusion 
length, 6 S , and influence factor, w) suffice to determine 
the geometrical features of a steady flow, such as the spa- 
tial distribution of mean velocity and diffusion, but two 
more are needed to introduce time dependence. The first 
is the mean rate of creating spots at the orifice (for sim- 
plicity, according to a Poisson process). In the DEM sim- 
ulation, particles exit a rate of mean rate of 4.40 x 10 3 s _1 , 
so spots carrying a typical volume V & 
be introduced at a mean rate of v s 
The second remaining spot parameter is the vertical drift 
speed, or, equivalently, the mean waiting time between 
spot displacements, At, which can be inferred from the 
drop in mean packing fraction during flow. In the DEM 
simulation, we find that there are initially 9, 400 parti- 
cles in the horizontal slice, 50d < z < 70<i, which drops 
to 8, 850 during flow. Choosing the spot waiting time to 
be At = 8.68 x 10 _4 s reproduces this decrease in den- 
sity in the spot simulation. The spot drift speed is thus 
v s = O.ld/ At = 115<i/s = 34.5cm/s, which is roughly ten 
times faster than typical particle speeds in Fig. 0] 



0.2051/p should 
2.15 x lOV 1 . 



IV. SPOT MODEL SIMULATION 

Having calibrated the five parameters (R s , b Sl w, v s , 
v s ), we can test the Spot Model by carrying out drainage 
simulations starting from the same static initial pack- 
ing as for the DEM simulations. For efficiency, a stan- 
dard cell method (also used in the parallel DEM code) is 
adapted for the spot simulations. The container is par- 
titioned into a grid of 10 x 3 x N z cells, each responsible 
for keeping track of the particles within it, with N z = 30 
for H = llOd and N z = 60 for H Q = 230d. When a spot 
moves, only the cells influenced by the spot need to be 
tested, and particles are transferred between cells when 
necessary. Without further optimization, the multiscale 
spot simulation runs over 100 times faster than the DEM 
simulation. 

The flow is initiated as spots are introduced uniformly 
at random positions on the orifice (at least R s away from 
the edges) at random times according to a Poisson pro- 
cess of rate v s . (The waiting time is thus an exponen- 
tial random variable of mean v~ x .) Once in the con- 
tainer, spots also move at random times with a mean 
waiting time, At = v s /Az s . Spot displacements in the 
bulk are chosen randomly from four displacement vec- 
tors, Ar s = (±Ax s , 0, Az s ) : (0, ±Ax s , Az s ), with equal 
probability, so spots perform directed random walks on 
a body-centered cubic lattice (with lattice parameter 
2Az s = 0.2<i). We make this simple choice to accelerate 
the simulation because more complicated, continuously 
distributed and/or smaller spot displacements with the 
same drift and diffusivity give very similar results. Spot 
centers are constrained not to come within d of a bound- 
ary, and once a spot reaches the top of the packing, it is 
removed from the simulation. More realistic models for 
the orifice, walls, and free surface are left for future work; 
here we focus on flowing packings in the bulk. 

The particles in the simulation move passively in re- 
sponse to spot displacements without any lattice con- 
straints. Although the influence of a spot can take a very 
general form j23| , the most important aspect is its length 
scale, so here we choose the simplest possible model in 
Eq. JJJ, where the spot influences particles uniformly in a 
sphere of radius R s . As shown in Fig.^a), we center the 
spot influence on the midpoint of its step, which seems 
the most consistent with the concept of moving intersti- 
tial volume from the initial to the final spot position. To 
be precise, when a spot moves from r s to r s + Ar s , all 
particles less than R s away from r s + Ar s /2 are displaced 
by — wAr s . 

To preserve realistic packings, we carry out a simple 
elastic relaxation after each spot-induced block motion, 
as in Fig. EJb). All particles within a radius R s + 2d 
of the midpoint of the spot displacement exert a soft- 
core repulsion on each other, if they begin to overlap. 
Rather than relaxing to equilibrium or integrating New- 
ton's laws, however, we use the simplest possible algo- 
rithm: Each pair of particles separated by less than d 
moves apart with identical and opposite displacements, 
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FIG. 5: (Color online) Time evolution of the random packing (from left to right) in DEM (top) and the spot simulation 
(bottom), for the H — 230d system, starting from the same initial state. Each image is a vertical slice through the center of 
the silo near the orifice well below the free surface. 



{d — r)a, for some constant a > 1. Similarly, a parti- 
cle within d/2 of a wall moves away by a displacement, 
(| — r)a. Particle positions are updated simultaneously 
once all pairings are considered, but those within the 
shell, R s + d < r < R s + 2d, more than one diameter 
away from the initial block motion, are held fixed to pre- 
vent long-range disruptions. 

It turns out that, due to the cooperative nature of 
Spot Model, only extremely small relaxation is required 
to enforce packing constraints, mainly near spot edges 
where some shear occurs. Here, we choose a = 0.8 and 
find that the displacements due to relaxation are typ- 



ically less than 25% of the initial block displacement, 
which is at the scale of 1/10,000 of a particle diameter: 
0.25wAr s %2x 10 -4 <i. Due to this tiny scale, the de- 
tails of the relaxation do not seem to be very important; 
we have obtained almost indistinguishable results with 
a — 0.6 and a = 1.0 and also with more complicated en- 
ergy minimization schemes. As such, we do not view the 
soft-core repulsion as introducing any new parameters. 
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V. RESULTS 

The spot and DEM simulations are compared using 
snapshots of all particle positions taken every 2 x 10 4 
time steps. As shown in Figure the agreement between 
the two simulations is remarkably good, considering the 
small number of parameters and physical assumptions 
in the Spot Model. It is clear a posteriori that the re- 
laxation step, in spite of causing only minuscule extra 
displacements, manages to produce reasonable packings 
during flow, while preserving the realistic description of 
the mean velocity and diffusion in the basic Spot Model. 
Only one parameter, b Sl is fitted to the mean flow, but 
we find that the entire velocity profile is accurately re- 
produced in the lower part of the container, as shown in 
Fig. 21 although the flow becomes somewhat more plug- 
like in DEM simulation higher in the container. Similarly, 
we fit w to the particle diffusion length in middle of the 
DEM simulation, b p = 2.86 x 10 -3 <i, without accounting 
for the elastic relaxation step, so it is reassuring that the 
same measurement in the spot simulation yields a similar 
value, b p = 2.73 x 10~ 3 d. 

The most surprising findings concern the agreement 
between the DEM and spot simulations for various mi- 
croscopic statistical quantities. First, we consider the ra- 
dial distribution function, g(r), which is the distribution 
of inter-particle separations, scaled to the same quantity 
in a ideal gas at the same density. For dense sphere pack- 
ings, the distribution begins with a large peak near r = d 
for particles in contact and smoothly connects smaller 
peaks at typical separations of more distant neighbors, 
while decaying to unity. As shown in Fig.UJa), the func- 
tions g(r) from the spot and DEM simulations are nearly 
indistinguishable, across the entire range of neighbors 
for the H Q = HOd system. This cannot be attributed 
entirely to the initial packing because each simulation 
evolves independently through substantial drainage and 
shearing. 

Next, we consider the three-body correlation function, 
#3(#), which gives the probability distribution for "bond 
angles" subtended by separation vectors to first neigh- 
bors (defined by separations less than the first minimum 
of g(r) at 1.38d). For sphere packings, gs(6) has a sharp 
peak at 60° for close-packed triangles, and another broad 
peak around 110 — 120° for larger crystal-like configura- 
tions. In Fig. EJb), we reach the same conclusion for 
#3(0) as for g(r): The spot and DEM simulations evolve 
independently from the initial packing to nearly indistin- 
guishable steady states. 

The striking agreement between the spot and DEM 
simulations seems to apply not only to structural, but 
also to dynamical, statistical quantities. Returning to 
Fig. 03 we see that the two simulations have very simi- 
lar spatial velocity correlations. Of course, the spot size, 
R s , in the Spot Model (without relaxation) was fitted 
roughly to the scale of the correlations in the DEM sim- 
ulation, but the multiscale spot simulation also manages 
to reproduce most of the fine structure of the correlation 
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FIG. 6: (Color online) Comparison of radial distribution func- 
tions (top) and bond angles (bottom) for particles in the re- 
gion — 15d < x < 15d, 15d < z < 45d for H Q = 110c? system. 
Three curves are shown on each graph, the first calculated 
from the initial static packing (common between the two sim- 
ulations), and the second and third calculated for over the 
range 1.04s < t < 1.40s. 

function. 

At much longer times, however, the random packings 
are no longer indistinguishable, as a small tendency for 
local close-packed ordering appears the spot simulation. 
As shown in Fig. [7[ the spot simulation develops en- 
hanced crystal-like peaks in g(r) at r = y/3d, 2d, 
The number of particles involved, however, is very small 
(~ 2%), and the effect seems to saturate, with no sig- 
nificant change between 8s and 16s. This is consistent 
with even longer spot simulations in systems with pe- 
riodic boundary conditions, which reach a similar, re- 
producible steady state (at the same volume fraction) 
from a variety of initial conditions 27]. In all cases, the 
spot algorithm never breaks down (e.g. due to jamming 
or instability), and unrealistic packings with overlapping 
particles are never created. 

The structure of the flowing steady state is fairly in- 
sensitive to various details of the spot algorithm. For ex- 
ample, changing the relaxation parameter (in the range 
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FIG. 7: (Color online) Evolution of the radial distribution 
function g(r) for H Q = 230d in the region — 15d < x < 15d, 
15d < z < 45gL The spot simulation (dashed curves) reaches 
a somewhat different steady state from the DEM simulation 
(solid curve), after a large amount of drainage has taken place. 

0.6 < a < 1.0), rescaling the spot size (by ±25%), and 
using a persistent random walk (for smoother spot tra- 
jectories), all have no appreciable effect on g(r). On the 
other hand, decreasing the vertical spot step size (in the 
range 0.025<i < Az < O.ld) tends to inhibit spurious lo- 
cal ordering and reduce the difference in g(r) between 
the spot and DEM simulations (e.g. measured by the 
norm). Therefore, our spot algorithm appears to "con- 
verge" with decreasing time step (and increasing com- 
putational cost), analogous to a finite- difference method, 
although this merits further study. 

VI. CONCLUSIONS 

Our results suggest that flowing dense random pack- 
ings have some universal geometrical features. This 
would be in contrast to static dense random packings, 
which suffer from ambiguities related to the degree of 
randomness and definitions of jamming @, 0- The simi- 



lar packing dynamics in spot and DEM simulations sug- 
gest that geometrical constraints dominate over mechan- 
ical forces in determining structural rearrangements, at 
least in granular drainage. Some form of the Spot Model 
may also apply to other granular flows and perhaps even 
to glassy relaxation, where localized, cooperative motion 
also occurs 

The Spot Model provides a simple framework for the 
multiscale modeling of liquids and glasses, analogous 
to dislocation dynamics in crystals. Our algorithm, 
which combines an efficient, "coarse-grained" simulation 
of spots with limited, local relaxation of particles, runs 
over 100 times faster than fully particle-based DEM for 
granular drainage. On current computers, this means 
that simulating one cycle of pebble-bed reactor [B| can 
take hours instead of weeks Q, although a general the- 
ory of spot motion in different geometries is still lacking. 
This may come from a stochastic formulation of Mohr- 
Coulomb plasticity, where spots perform random walks 
along slip lines of incipient failure [2^, which could, in 
principle, be applied to different materials by changing 
the yield criterion. Alternatively, a multiscale model for 
supercooled molecular liquids could involve spots mov- 
ing along chains of dynamic facilitation [2(J, u$ • In any 
case, we have demonstrated that dense random-packing 
dynamics can be driven entirely by the motion of simple, 
collective excitations. 
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